//===================================
// G.Tatishvili
// PNNL
// hb-analysis code. Version 1.0.0
//===================================

#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>

#include "CLHEP/Vector/LorentzVector.h"
#include "helix/Helix.h"

#include "event/BelleEvent.h"
#include "tuple/BelleTupleManager.h"

#include "basf/module.h"
#include "basf/module_descr.h"

#include "panther/panther.h"
#include BELLETDF_H
#include HEPEVT_H
#include MDST_H
#include MDST_OBS_H
#include EVTCLS_H
#include EVTVTX_H
#include LEVEL4_H
#include TRG_H

#include "ip/IpProfile.h"
#include "kid/atc_pid.h"
#include "eid/eid.h"
#include "mdst/mdst.h"
#include "mdst/Muid_mdst.h"
#include "mdst/findKs.h"

#include "kfitter/kvertexfitter.h"
#include "kfitter/kmassfitter.h"
#include "kfitter/kmassvertexfitter.h"
#include "kfitter/kmakemother.h"

#include "particle/Particle.h"
#include "particle/utility.h"
#include "particle/gammac.h"

#include "benergy/BeamEnergy.h"

#if defined(BELLE_NAMESPACE)
namespace Belle {
#endif

extern "C" {
  float prob_(const float&, const int&);
}


#define E_MASS   0.000511
#define MU_MASS  0.105658
#define PI_MASS  0.139570
#define K_MASS   0.493677
#define P_MASS   0.938272

#define LGNDR2(X) ((3.*X*X-1.)/2.)
#define LGNDR3(X) ((5.*X*X*X-3.*X)/2.)
#define LGNDR4(X) ((35.*X*X*X*X-30.*X*X+3.)/8.)
#define LGNDR5(X) ((63.*X*X*X*X*X-70.*X*X*X+15.*X)/8.)
#define LGNDR6(X) ((231.*X*X*X*X*X*X-315.*X*X*X*X+105.*X*X-5.)/16.)
#define LGNDR7(X) ((429.*X*X*X*X*X*X*X-693.*X*X*X*X*X+315.*X*X*X-35.*X)/16.)
#define LGNDR8(X) ((6435.*X*X*X*X*X*X*X*X-12012.*X*X*X*X*X*X+6930.*X*X*X*X-1260.*X*X+35)/128.)

  //fix below OR remove boost from HepParticles
  /*
  double E_HER = BeamEnergy::E_HER();
  double E_LER = BeamEnergy::E_LER();
  double cross_angle = BeamEnergy::Cross_angle();

  HepLorentzVector p4_her(sqrt(E_HER*E_HER-E_MASS*E_MASS)*sin(cross_angle),
                          0.0,
                          +sqrt(E_HER*E_HER-E_MASS*E_MASS)*cos(cross_angle),
                          E_HER);
  //p4_her.rotateY((float) cross_angle);
  HepLorentzVector p4_ler(0.0,
                          0.0,
                          -sqrt(E_LER*E_LER-E_MASS*E_MASS),
                          E_LER);
  HepLorentzVector P_BEAM = (p4_her + p4_ler);
  
  */

  HepPoint3D   IP(0,0,0);
  HepSymMatrix IPerr(3,0);

  atc_pid kid(3,1,5,3,2);
  atc_pid prid(3,1,5,4,2);




  class HepParticle {
  public:

    HepParticle() {charged=NULL;gamma=NULL;pi0=NULL;vee2=NULL;mode=-1;charge=0;
      child[0]=child[1]=child[2]=child[3]=NULL;}

    HepParticle(Mdst_charged& trk, double mass, int type=-1) { // Constructer for charged particles
      if (type>=0 && type<=4) 
	{
	  Mdst_trk_fit& trkFit = trk.trk().mhyp(type);
	  
	  HepVector a(5,0);
	  a[0] = trkFit.helix(0);
	  a[1] = trkFit.helix(1);
	  a[2] = trkFit.helix(2);
	  a[3] = trkFit.helix(3);
	  a[4] = trkFit.helix(4);
	  
	  HepPoint3D pivot(trkFit.pivot(0), trkFit.pivot(1), trkFit.pivot(2));
	  Helix ltrk(pivot, a);
	  ltrk.pivot(HepPoint3D(0.,0.,0.));

	  P4 = ltrk.momentum(0., mass);
	}
      else 
	{
	  HepLorentzVector P ( trk.px(), trk.py(), trk.pz(),
		      sqrt ( trk.px()*trk.px() + trk.py()*trk.py() +
			     trk.pz()*trk.pz() + mass*mass ) );
	  P4    = P;
	}

      gamma=NULL;pi0=NULL;vee2=NULL;
      child[0]=child[1]=child[2]=child[3]=NULL;
      charged = &trk;
      charge = (int)trk.charge();
      mode=-1;

      //            P4_boost = P4;
      //    P4_boost.boost(-P_BEAM.boostVector());
    }

    HepParticle(Mdst_gamma& gam) { // Constructer for gamma
      HepLorentzVector P ( gam.px(), gam.py(), gam.pz(),
		  sqrt ( gam.px()*gam.px() + gam.py()*gam.py() +
			 gam.pz()*gam.pz()) );
      P4    = P;
      charged=NULL;pi0=NULL;vee2=NULL;
      child[0]=child[1]=child[2]=child[3]=NULL;
      mode=-1;
      gamma = &gam;
      charge = 0;

      //      P4_boost = P4;
      //P4_boost.boost(-P_BEAM.boostVector());

      //      cout << "P_BEAM 1 " << P_BEAM << endl;

    }

    HepParticle(Mdst_pi0& mdstpi0) { // Constructer for pi0
      HepLorentzVector P ( mdstpi0.px(), mdstpi0.py(), mdstpi0.pz(),
		  mdstpi0.energy() );
      P4  = P;
      charged=NULL;gamma=NULL;vee2=NULL;
      child[0]=child[1]=child[2]=child[3]=NULL;
      mode=0;
      pi0 = &mdstpi0;
      charge = 0;

      //      P4_boost = P4;
      //P4_boost.boost(-P_BEAM.boostVector());
    }
    HepParticle(Mdst_vee2& mdstvee2) { // Constructer for Ks
      HepLorentzVector P ( mdstvee2.px(), mdstvee2.py(), mdstvee2.pz(),
		  mdstvee2.energy() );
      P4  = P;
      charged=NULL;gamma=NULL;pi0=NULL;
      child[0]=child[1]=child[2]=child[3]=NULL;
      mode=0;
      vee2 = &mdstvee2;
      charge = 0;

      //      P4_boost = P4;
      //P4_boost.boost(-P_BEAM.boostVector());
    }
    HepParticle(HepParticle& ch1, HepParticle& ch2) { // 2-body Constructer
      P4=ch1.P4+ch2.P4;
      charge=ch1.charge+ch2.charge;
      child[0]=&ch1;
      child[1]=&ch2;
      child[2]=NULL;
      child[3]=NULL;
      charged=NULL;gamma=NULL;pi0=NULL;vee2=NULL;
      mode=0;


      //      P4_boost = P4;
      //P4_boost.boost(-P_BEAM.boostVector());
    }
    HepParticle(HepParticle& ch1, HepParticle& ch2,
		HepParticle& ch3) {                   // 3-body Constructer
      P4=ch1.P4+ch2.P4+ch3.P4;
      charge=ch1.charge+ch2.charge+ch3.charge;
      child[0]=&ch1;
      child[1]=&ch2;
      child[2]=&ch3;
      child[3]=NULL;
      charged=NULL;gamma=NULL;pi0=NULL;vee2=NULL;
      mode=0;

      //      P4_boost = P4;
      //P4_boost.boost(-P_BEAM.boostVector());
    }
    HepParticle(HepParticle& ch1, HepParticle& ch2,
		HepParticle& ch3, HepParticle& ch4) { // 4-body Constructer
      P4=ch1.P4+ch2.P4+ch3.P4+ch4.P4;
      charge=ch1.charge+ch2.charge+ch3.charge+ch4.charge;
      child[0]=&ch1;
      child[1]=&ch2;
      child[2]=&ch3;
      child[3]=&ch4;
      charged=NULL;gamma=NULL;pi0=NULL;vee2=NULL;
      mode=0;

      //      P4_boost = P4;
      //P4_boost.boost(-P_BEAM.boostVector());
    }

  ~HepParticle() {};

    HepParticle & operator = (const HepParticle &p) { // Assigment
      P4 = p.P4;
      P4_boost = p.P4_boost;
      child[0] = p.child[0];
      child[1] = p.child[1];
      child[2] = p.child[2];
      child[3] = p.child[3];
      charged = p.charged;
      gamma = p.gamma;
      pi0 = p.pi0;
      vee2 = p.vee2;
      vertex = p.vertex;
      charge = p.charge;
      mode = p.mode;
    }
    


    void GetImpactParameters(double *dr, double *dz, int *rphi, int *z, int *ax, int *st, int t, HepPoint3D IPfrom = IP) {

      *dr = 1000.0;
      *dz = 1000.0;
      *rphi = 0;
      *z = 0;
      *ax = 0;
      *st = 0;

      if(charged->trk()){
	Mdst_trk_fit& trkFit = charged->trk().mhyp(t);
	//	cout << "SVD hitsrp " << trkFit.nhits(3) << endl;
	//cout << "SVD hitsz " << trkFit.nhits(4) << endl;
	if(trkFit){
	  HepVector a(5,0);
	  a[0] = trkFit.helix(0);
	  a[1] = trkFit.helix(1);
	  a[2] = trkFit.helix(2);
	  a[3] = trkFit.helix(3);
	  a[4] = trkFit.helix(4);

	  HepPoint3D pivot(trkFit.pivot(0), trkFit.pivot(1), trkFit.pivot(2));
	  Helix ltrk(pivot, a);
	  ltrk.pivot(IPfrom);

	  *dr = ltrk.dr();
	  *dz = ltrk.dz();
	  *rphi = trkFit.nhits(3);
	  *z = trkFit.nhits(4);
	  *ax = trkFit.nhits(0);
	  *st = trkFit.nhits(1);
	}
      }
    }
     
  // General informations
  HepLorentzVector P4;    		    	    	                // Four momentum in Lab-frame 
  HepLorentzVector P4_boost;		    	    	                // Four momentum in CM-frame
  HepParticle *child[4];                                        // Children
  Mdst_charged *charged;	    	    	                // Mdst banks for Pi,K
  Mdst_gamma *gamma;		    	    	                // Mdst banks for Gamma
  Mdst_pi0 *pi0;		    	    	                // Mdst banks for pi0
  Mdst_vee2 *vee2;  	    	    	    	                // Mdst banks for Vee2
  HepPoint3D vertex;                                            // vertex point
  int mode; 	    	    	    	    	                // mode code
  int charge;                                                   // charge of the particle
  };

  // Class definition
  class hb : public Module {
  public:
    hb ( void ){ eventcount=savecount=0;};
    ~hb ( void ){};
    void init ( int * ){};
    void term ( void ){
      printf("Final - %d/%d\n",savecount,eventcount); 
    };
    void disp_stat ( const char* ) {} ;
    void hist_def ( void );
    void event ( BelleEvent*, int* );
    void begin_run ( BelleEvent*, int * ){    
      eid::init_data();

      //get Beam Energy info
      BeamEnergy::begin_run();


      // get IP profile data from database
      IpProfile::begin_run();

      // dump IP profile data to STDOUT
      IpProfile::dump();    
    };
    void end_run ( BelleEvent*, int * ){};
    void other ( int*, BelleEvent*, int* ){}; 

private:
    BelleTuple *m_tpl, *m_tplpp;
    BelleTuple *m_new;
    BelleHistogram *h_hbmass,*h_hbprimemass;
    BelleHistogram *h_etabmass,*h_etabprimemass;
    BelleHistogram *h_Y5Smass,*h_Y4Smass,*h_Y3Smass,*h_Y2Smass,*h_Y1Smass;
    BelleHistogram *h_vtx, *h_vtxsv, *h_vtxprob;
    BelleHistogram *h_ntrk, *h_ntrksv, *h_dimu;

    void fwjet2(int NTrk, float *PTrk, float EBeam, 
		float *H0, float *R1, float *R2, float *R3, float *R4);


   
  int MCTruthInfo(HepParticle& Rec);
  int eventcount,savecount;
};

  //  extern float PROB(float chisq, int dof);

  extern "C" Module_descr *mdcl_hb ()
  {
    hb *module = new hb;
    Module_descr *dscr = new Module_descr ( "hb", module );
    IpProfile::define_global(dscr);
    return dscr;
  }
                                          
// Member Functions

// hist_def function
  void hb::hist_def ( void )
  {  
    extern BelleTupleManager *BASF_Histogram;
    BelleTupleManager& tm = *BASF_Histogram;
    BelleTupleManager& tmpp = *BASF_Histogram;
    BelleTupleManager& tmnew = *BASF_Histogram;
    
    /* Large nt for MC studies

    m_tpl = tm.ntuple ("Event", "ntr npi0 ngam npostr nnegtr qpi1 qpi2 \
drpi1 drpi2 dzpi1 dzpi2 kprobpi1 kprobpi2  \
pxpi1 pypi1 pzpi1 plabpi1 pxpi2 pypi2 pzpi2 plabpi2 pxg pyg \
pzg eg egboost recm2pi recmg2pi m2pi angle2pi ebeam expno runno evtno \
ncls nclsb psum esum evis evtpz hjmass thrust evtr2");

    m_tplpp = tmpp.ntuple ("Event pipi", "ntr npi0 ngam npostr nnegtr qpi1 qpi2 \
drpi1 drpi2 dzpi1 dzpi2 kprobpi1 kprobpi2 \
pxpi1 pypi1 pzpi1 plabpi1 pxpi2 pypi2 pzpi2 plabpi2 \
recm2pi m2pi angle2pi ebeam expno runno evtno \
ncls nclsb psum esum evis evtpz hjmass thrust evtr2");
*/

    //small ntupes for data
    m_tpl = tm.ntuple ("Event", "ntr npi0 ngam npostr nnegtr qpi1 qpi2 \
rp1 rz1 rp2 rz2 drpi1 drpi2 dzpi1 dzpi2 kprobpi1 kprobpi2 plabpi1 plabpi2 \
pzg eg egboost recm2pi recmg2pi m2pi angle2pi ebeam expno runno evtno \
ncls nclsb psum esum evis evtpz hjmass thrust evtr2 \
chisq chisqsv chisq2pi");

    m_tplpp = tmpp.ntuple ("Event pipi", "ntr npi0 ngam npostr nnegtr qpi1 qpi2 \
rp1 rz1 rp2 rz2 drpi1 drpi2 dzpi1 dzpi2 kprobpi1 kprobpi2 plabpi1 plabpi2 \
recm2pi m2pi angle2pi ebeam expno runno evtno \
ncls nclsb psum esum evis evtpz hjmass thrust evtr2 \
chisq chisqsv chisq2pi");

    /*        m_new = tmnew.ntuple ("Event new", "ntr npi0 ngam0 \
expno runno evtno energy prob2pi probev \
svd2pi angle2pi kprobpi1 kprobpi2 pprobpi1 pprobpi2 eid1 eid2 muid1 muid2 evtr2 recm2pi m2pi m2mu recmg2pi \
angleg eg egboost pqpi1 pqpi2 anglelabg dr1 dr2 dz1 dz2 rp1 rp2 rz1 rz2 ax1 ax2 st1 st2 recmpi1 recmpi2 npostr nnegtr mchrg npipi ngam1 ngam2 ngam3 ngam4 ngam5 ngam6 ngam7 ngam8 ngam9");
    */
    // slightly reduced m_new below - removed svd2pi, ax1, ax2, st1, st2, probev, npostr, nnegtr, mchrg
        m_new = tmnew.ntuple ("Event new", "ntr npi0 ngam0 \
expno runno evtno energy prob2pi \
angle2pi kprobpi1 kprobpi2 pprobpi1 pprobpi2 eid1 eid2 muid1 muid2 evtr2 recm2pi m2pi m2mu recmg2pi \
angleg eg egboost pqpi1 pqpi2 anglelabg dr1 dr2 dz1 dz2 rp1 rp2 rz1 rz2 recmpi1 recmpi2 npipi ngam1 ngam2 ngam3 ngam4 ngam5 ngam6 ngam7 ngam8 ngam9");

    h_hbmass  = BASF_Histogram->histogram("hb MC mass", 100, 9.0, 10.0);
    h_etabmass  = BASF_Histogram->histogram("etab MC mass", 100,9.0, 10.0);
    h_hbprimemass  = BASF_Histogram->histogram("hb' MC mass", 100, 9.0, 10.0);
    h_etabprimemass  = BASF_Histogram->histogram("etab' MC mass", 100,9.0, 10.0);
    h_Y5Smass = BASF_Histogram->histogram("Y(5S) mass", 100, 9.0, 11.0);
    h_Y4Smass = BASF_Histogram->histogram("Y(4S) mass", 100, 9.0, 11.0);
    h_Y3Smass = BASF_Histogram->histogram("Y(3S) mass", 100, 9.0, 11.0);
    h_Y2Smass = BASF_Histogram->histogram("Y(2S) mass", 100, 9.0, 11.0);
    h_Y1Smass = BASF_Histogram->histogram("Y(1S) mass", 100, 9.0, 11.0);

    h_vtx = BASF_Histogram->histogram("vtx chisq", 102, -2, 100);
    h_vtxsv = BASF_Histogram->histogram("vtx chisq with sv", 102, -2, 100);
    h_ntrk = BASF_Histogram->histogram("vtx trk", 21, -0.5, 20.5);
    h_ntrksv = BASF_Histogram->histogram("vtx trk with sv", 21, -0.5, 20.5);
    h_vtxprob = BASF_Histogram->histogram("vtx prob with sv", 100, 0., 1.);

    h_dimu = BASF_Histogram->histogram("DiMuon Mass", 100, 9., 11.);

  }

// event function
void hb::event ( BelleEvent* evptr, int* status )
{
  // Get managers for MDST banks
  Mdst_charged_Manager& ChgMgr = Mdst_charged_Manager::get_manager();
  Mdst_gamma_Manager& GamMgr = Mdst_gamma_Manager::get_manager();
  Mdst_pi0_Manager& Pi0Mgr = Mdst_pi0_Manager::get_manager();
  Mdst_vee2_Manager& Vee2Mgr = Mdst_vee2_Manager::get_manager();
  Belle_event_Manager& EvtMgr = Belle_event_Manager::get_manager();
  Gen_hepevt_Manager& HepevtMgr = Gen_hepevt_Manager::get_manager();  
  Mdst_ecl_Manager& EclMgr = Mdst_ecl_Manager::get_manager();
  Mdst_ecl_aux_Manager& EclAuxMgr = Mdst_ecl_aux_Manager::get_manager();
  Evtcls_hadronic_flag_Manager& HadFlagMgr = Evtcls_hadronic_flag_Manager::get_manager();
  Evtcls_flag_Manager& FlagMgr = Evtcls_flag_Manager::get_manager();
  Evtcls_hadron_info_Manager& HadInfoMgr = Evtcls_hadron_info_Manager::get_manager();
  L4_summary_Manager& L4SumMgr = L4_summary_Manager::get_manager();
  Evtvtx_primary_vertex_Manager& VtxMgr = Evtvtx_primary_vertex_Manager::get_manager();
  Rectrg_summary3_Manager& TrgSumMgr = Rectrg_summary3_Manager::get_manager();

  *status=0;


  // set IP and error
  int IPUsable = 0;
  if(IpProfile::usable()){
    IP    = IpProfile::position(1);
    IPerr = IpProfile::position_err_b_life_smeared(1);
    IPUsable = 1;
  }else {
    IP    = HepPoint3D(0,0,0);
    IPerr = HepSymMatrix(3,0);
  }


  //--------------------------------

  /*  float e_cms = 10.8667;

  double E_HER_0 = 7.998213;
  double E_LER_0 = 3.499218;
  double E_BEAM  = e_cms/2;
  double E_HER = (E_BEAM/5.29)*E_HER_0;
  double E_LER = (E_BEAM/5.29)*E_LER_0;
  double cross_angle = 0.022;
  double MASS_E = 0.000511;
  */

  double E_BEAM = BeamEnergy::E_beam_corr();
  double E_HER = BeamEnergy::E_HER();
  double E_LER = BeamEnergy::E_LER();
  double cross_angle = BeamEnergy::Cross_angle();

  HepLorentzVector p4_her(0.0,
                          0.0,
                          +sqrt(E_HER*E_HER-E_MASS*E_MASS),
                          E_HER);
  p4_her.rotateY(cross_angle);
  HepLorentzVector p4_ler(0.0,
                          0.0,
                          -sqrt(E_LER*E_LER-E_MASS*E_MASS),
                          E_LER);
  HepLorentzVector P_BEAM = (p4_her + p4_ler);
  
  //  cout << "P_BEAM 2 " << P_BEAM << endl;

  //--------------------------------


  int Necl_barrel = 0;
  for(std::vector<Mdst_ecl>::iterator it = EclMgr.begin();
      it != EclMgr.end(); it++){
    Mdst_ecl& cls = *it;

    if(cls.quality() == 0 && cls.energy() > 0.1 && //E_min_cut 
       cos(cls.theta()) > -0.7 && // ECLcos_min
       cos(cls.theta()) <  0.9 )  // ECLcos_max
      Necl_barrel++;
  }


  int evt_Ntrk = -1, evt_Ncls = -1;
  double evt_Psum, evt_Esum, evt_Evis, evt_Pz, evt_HeavyJetMass, evt_Thrust, evt_R2;
  evt_Psum = evt_Esum = evt_Evis = evt_Pz = evt_HeavyJetMass = evt_Thrust = evt_R2 = -1.;
  
  if (HadInfoMgr.count()>0) {
    evt_Ntrk = HadInfoMgr[0].Ntrk(); 
    evt_Ncls = HadInfoMgr[0].Ncls();
    evt_Psum = HadInfoMgr[0].Psum();
    evt_Esum = HadInfoMgr[0].Esum();
    evt_Evis = HadInfoMgr[0].Evis();
    evt_Pz   = HadInfoMgr[0].Pz();
    evt_HeavyJetMass = HadInfoMgr[0].HeavyJetMass();
    evt_Thrust    = HadInfoMgr[0].Thrust();
    evt_R2    = HadInfoMgr[0].R2();
  }

// Define Track list
  std::vector<HepParticle> TrkList,GamList,Pi0List, MupList, MumList;
  std::vector<Particle> PartTrkList;

  int Npos = 0;
  int Nneg = 0;

  HepLorentzVector P_CHG;

// Loop over charged tracks
  for(std::vector<Mdst_charged>::iterator it = ChgMgr.begin(); it != ChgMgr.end(); ++it)
    {
      Mdst_charged& Charged = *it;
      Particle Partpl(Charged,"PI+"); 
      Particle Partmi(Charged,"PI-"); 


      if( !good_charged(*it) ) continue;

      if (Charged.charge()>0)
	{ Npos++;            
	}
      if (Charged.charge()<0) 
	{Nneg++;     
	}

      HepParticle part(Charged,PI_MASS,2); 

      Muid_mdst muID(Charged);
	

     
      

      Mdst_trk_fit &trk_fit = Charged.trk().mhyp(2);

      //            if (Charged.trk().quality() == 0
	  //	  	  && trk_fit.nhits(3)>0
	  //&& trk_fit.nhits(4)>1
      //      	  )
	{
	  //	  cout << "PI_MASS="<<PI_MASS<<endl;
	  TrkList.push_back(part);
	  if (Charged.charge()>0) {PartTrkList.push_back(Partpl);}
	  if (Charged.charge()<0) {PartTrkList.push_back(Partmi);}

	  P_CHG = P_CHG + part.P4;

	  if (muID.Chi_2()>0){
	    double mu_id = muID.Muon_likelihood();
	    //	    cout << "muid " << mu_id << endl;
	    	    if (mu_id > 0.9)
	      {
		if (Charged.charge()>0) {MupList.push_back(part);}
		if (Charged.charge()<0) {MumList.push_back(part);}
	      }
	  }
	}
    }



// Loop over gamma 
  for(std::vector<Mdst_gamma>::iterator it = GamMgr.begin(); it != GamMgr.end(); ++it)
    {    
      Mdst_gamma& Gamma = *it;           
      HepParticle part(Gamma);
      double cos_theta = part.P4.pz()/part.P4.rho();
      Mdst_ecl_aux &aux = EclAuxMgr( Panther_ID(Gamma.ecl().get_ID()) );
      double e9oe25 = aux.e9oe25();
	
      //      cout << "e9oe25=" << e9oe25 << endl;

      if ((Gamma.ecl().match() == 0 || Gamma.ecl().match() == 2)
	  //&&		 sqrt ( Gamma.px()*Gamma.px() +Gamma.py()*Gamma.py() +
	  //	   Gamma.pz()*Gamma.pz())> 0.5)
	  )
      {    
	if (e9oe25 < 0.95) continue;
	if (cos_theta > 0.8501 && part.P4.e() < 0.100) continue;
	if (-0.6643 < cos_theta && cos_theta < 0.8501 && part.P4.e() < 0.050) continue;
	if (-0.6643 < cos_theta && cos_theta < 0.8501 && part.P4.e() < 0.1) continue;
	if (cos_theta < -0.6643 && part.P4.e() < 0.100) continue;


	//		if (cos_theta > 0.8501) continue; //not forward barrel


	HepParticle part(Gamma);
	GamList.push_back(part);
	}
    }


  int ntrk=TrkList.size()+GamList.size(),trkcount=0;
  float ptrk[ntrk*3],H0,R1,R2,R3,R4;

  for(std::vector<HepParticle>::iterator it = TrkList.begin();
      it != TrkList.end(); ++it){

    HepParticle& Charged = *it;

    ptrk[trkcount*3]=Charged.P4_boost.px();
    ptrk[trkcount*3+1]=Charged.P4_boost.py();
    ptrk[trkcount*3+2]=Charged.P4_boost.pz();
    trkcount++;
  }

  for(std::vector<HepParticle>::iterator it = GamList.begin();
      it != GamList.end(); ++it){

    HepParticle& Gamma = *it;

    ptrk[trkcount*3]=Gamma.P4_boost.px();
    ptrk[trkcount*3+1]=Gamma.P4_boost.py();
    ptrk[trkcount*3+2]=Gamma.P4_boost.pz();
    trkcount++;
  }

  fwjet2(ntrk,ptrk,E_BEAM,&H0,&R1,&R2,&R3,&R4);

  int pi0_count = 0;

  for(std::vector<Mdst_pi0>::iterator it = Pi0Mgr.begin();
      it != Pi0Mgr.end(); ++it){        
  
    Mdst_pi0& Pi0 = *it;    

    if (fabs(0.134-Pi0.mass())>=0.0054*3) continue;

    HepParticle part(Pi0);  // Create a Pi0 as HepParticle  
    part.child[0] = part.child[1] = NULL;
    for(std::vector<HepParticle>::iterator it = GamList.begin();
	it != GamList.end(); ++it) {

      HepParticle& ga = *it;
      if (Pi0.gamma(0) == (*ga.gamma)) part.child[0] = &ga;
      if (Pi0.gamma(1) == (*ga.gamma)) part.child[1] = &ga;

    }
    if (part.child[0]!=NULL && part.child[1]!=NULL)
      {
	Pi0List.push_back(part);  
	pi0_count++;
      }
  }

  //make gamma list that is not included among pi0 photons

  std::vector<HepParticle> myGamList;


  for(std::vector<HepParticle>::iterator it0 = GamList.begin();
      it0 != GamList.end(); ++it0){
    
    HepParticle& Gamma0 = *it0;
    int mypi0_count = 0;

    
    for(std::vector<HepParticle>::iterator it1 = GamList.begin();
	it1 != GamList.end(); ++it1){
      
      HepParticle& Gamma1 = *it1;

      HepLorentzVector P_Pi0 = Gamma0.P4+Gamma1.P4;
      
      double Mass_Pi0 = P_Pi0.m();
      
      if (fabs(0.134-Mass_Pi0)<=0.0054*3){mypi0_count++;}
      
    }

     if (mypi0_count==0)
      {
	HepParticle part(Gamma0);
	myGamList.push_back(part);
	/*if (Gamma0.P4.e()>=0.100&&Gamma0.P4.e()<0.150){
	  if (my_ngam[0]<9){ my_ngam[0]=my_ngam[0]+1;}}
	if (Gamma0.P4.e()>=0.150&&Gamma0.P4.e()<0.200){
	  if (my_ngam[1]<9){ my_ngam[1]=my_ngam[1]+1;}}
	if (Gamma0.P4.e()>=0.200&&Gamma0.P4.e()<0.250){
	  if (my_ngam[2]<9){ my_ngam[2]=my_ngam[2]+1;}}
	if (Gamma0.P4.e()>=0.250&&Gamma0.P4.e()<0.300){
	  if (my_ngam[3]<9){ my_ngam[3]=my_ngam[3]+1;}}
	if (Gamma0.P4.e()>=0.300&&Gamma0.P4.e()<0.350){
	  if (my_ngam[4]<9){ my_ngam[4]=my_ngam[4]+1;}}
	if (Gamma0.P4.e()>=0.350&&Gamma0.P4.e()<0.400){
	  if (my_ngam[5]<9){ my_ngam[5]=my_ngam[5]+1;}}
	if (Gamma0.P4.e()>=0.400&&Gamma0.P4.e()<0.450){
	  if (my_ngam[6]<9){ my_ngam[6]=my_ngam[6]+1;}}
	if (Gamma0.P4.e()>=0.450&&Gamma0.P4.e()<0.500){
	  if (my_ngam[7]<9){ my_ngam[7]=my_ngam[7]+1;}}
	if (Gamma0.P4.e()>=0.500&&Gamma0.P4.e()<0.550){
	  if (my_ngam[8]<9){ my_ngam[8]=my_ngam[8]+1;}}
	if (Gamma0.P4.e()>=0.550){
	  if (my_ngam[9]<9){ my_ngam[9]=my_ngam[9]+1;}}*/
      }
  }
  //end attempt to make new gamma list

  
    

  Gen_hepevt_Manager &hepmgr = Gen_hepevt_Manager::get_manager();
  for (std::vector<Gen_hepevt>::iterator it = hepmgr.begin();
       it != hepmgr.end(); ++it) {
    
    if (abs(it->idhep()) == 10553){
      h_hbmass->accumulate ( (float) it->M(), 1.0);
    }
    if (abs(it->idhep()) == 551){
      h_etabmass->accumulate ( (float) it->M(), 1.0);
    }
    if (abs(it->idhep()) == 110553){
      h_hbprimemass->accumulate ( (float) it->M(), 1.0);
    }
    if (abs(it->idhep()) == 100551){
      h_etabprimemass->accumulate ( (float) it->M(), 1.0);
    }
    if (abs(it->idhep()) == 553){
      h_Y1Smass->accumulate ( (float) it->M(), 1.0);
    }
    if (abs(it->idhep()) == 100553){
      h_Y2Smass->accumulate ( (float) it->M(), 1.0);
    }
    if (abs(it->idhep()) == 200553){
      h_Y3Smass->accumulate ( (float) it->M(), 1.0);
    }
    if (abs(it->idhep()) == 300553){
      h_Y4Smass->accumulate ( (float) it->M(), 1.0);
    }
    if (abs(it->idhep()) == 9000553){
      h_Y5Smass->accumulate ( (float) it->M(), 1.0);
    }
  }

  if (Npos <1 || Nneg <1 ) return;

  // vertex in this loop


  kvertexfitter kv, kvsv;
  int i=0, isv=0;
  std::vector<HepParticle>::iterator tit1 = TrkList.begin();
  for(std::vector<Particle>::iterator it1 = PartTrkList.begin(); it1 != PartTrkList.end()-1; ++it1)
    {

      Particle& Part= *it1;
      HepParticle& Pi1  = *tit1;
      double dr1,dz1;
      int nrp1, nz1;
      int nax1, nst1;
	  
      Pi1.GetImpactParameters(&dr1,&dz1,&nrp1,&nz1,&nax1,&nst1,0);      
      //	  Pi1.GetImpactParameters(&dr1,&dz1,&nrp1,&nz1,0);
      
      //if (enough_svd_hit(Part))
      {
	//	cout<< "dr,dz,nrpm,nz = " << dr1 << " " << dz1 << " " << nrp1 << " " << nz1 << endl;
	addTrack2fit(kv, Part);
	i++;
      }
      
      /*      Mdst_trk_fit &trk_fit = Part.mdstCharged().trk().mhyp(2);
      
      nrp1 = trk_fit.nhits(3);
      nz1 = trk_fit.nhits(4);*/

      //      cout << "probev crit " << dr1 << " " << dz1 << " " << nrp1 << " " << nz1 << " " << nax1 << " " << nst1;

      if (nrp1>2 && nz1 > 2 && (nax1+nst1)>19 && fabs(dr1)<0.03 && fabs(dz1)<1)
	{
	  //	  cout<< "dr,dz,nrpm,nz = " << dr1 << " " << dz1 << " " << nrp1 << " " << nz1 << endl;
	  addTrack2fit(kvsv, Part);
	  isv++;
	}
      
            tit1++;

    }

  if (IpProfile::usable()) addTube2fit(kv);
  if (IpProfile::usable()) addTube2fit(kvsv);
  
  int fitstatus_sv = -1;

  kv.fit();
  fitstatus_sv = kvsv.fit();
  /*    	      cout << "trk # " << i << endl;
	      cout << "kvsv.fit() "  << kvsv.fit() <<endl;
	      cout << "kv.vertex() "  << kv.vertex() <<endl;
	      cout << "kv.errVertex() "  << kv.errVertex() <<endl;
	      cout << "kv.chisq() "  << kv.chisq() <<endl;

	      cout << "trk sv# " << isv << endl;
	      cout << "kvsv.fit() "  << kvsv.fit() <<endl;
	      cout << "kvsv.vertex() "  << kvsv.vertex() <<endl;
	      cout << "kvsv.errVertex() "  << kvsv.errVertex() <<endl;
	      cout << "kvsv.chisq() "  << kvsv.chisq() <<endl;
  */  
	      float fchisq = kv.chisq();
	      float fchisqsv = -10;
	      if (isv>0) fchisqsv=kvsv.chisq();

	      //	      	      cout << "prob = " << prob_(fchisqsv, 2*isv) << endl;
	      h_vtx->accumulate ( fchisq, 1.0);
	      h_ntrk->accumulate ( (float) i, 1.0);
	      h_vtxsv->accumulate ( fchisqsv, 1.0);
	      h_ntrksv->accumulate ( (float) isv, 1.0);
	      float probev = -1;
	      if (isv>1 && fitstatus_sv==0) 
		{  if (IpProfile::usable()) probev = prob_(fchisqsv, 2*isv-1);
		  if (!IpProfile::usable()) probev = prob_(fchisqsv, 2*isv-3);
		  h_vtxprob->accumulate ( probev, 1.0);
		}

  //now for main loop

	      int my_npipi[10];
	      int pipi_bin = 0;
	      
	      for (pipi_bin=0;pipi_bin<10;++pipi_bin){
		my_npipi[pipi_bin]=0;
	      }


	      std::vector<Mdst_charged>::iterator eit1 = ChgMgr.begin();
	      std::vector<Particle>::iterator nit1 = PartTrkList.begin();
      for(std::vector<HepParticle>::iterator it1 = TrkList.begin();it1 != TrkList.end()-1; ++it1)
	{
	  std::vector<Mdst_charged>::iterator eit2 = eit1+1;
	  std::vector<Particle>::iterator nit2 = nit1+1;

	  for(std::vector<HepParticle>::iterator it2 = it1+1; it2 != TrkList.end(); 
	      ++it2)
	    {
	      

	      
	  HepParticle& Pi1  = *it1;
	  HepParticle& Pi2  = *it2;

	  Particle& Part1= *nit1;
	  Particle& Part2= *nit2;

	  Mdst_charged& el1 = *eit1;
	  Mdst_charged& el2 = *eit2;

	  //	  Mdst_charged& trk1 = *it1;
	  //Mdst_charged& trk2 = *it2;
	    
	  //Mdst_trk_fit &trk_fit1 = trk1.mdstCharged().trk().mhyp(2);
	  //Mdst_trk_fit &trk_fit2 = trk2.mdstCharged().trk().mhyp(2);

	  double dr1,dz1,dr2,dz2;
	  int nrp1, nz1, nrp2, nz2, nax1, nst1, nax2, nst2;

	  Muid_mdst muID1(el1), muID2(el2);
	  double mu_id1 = -1;
	  double mu_id2 = -1;
	  if (muID1.Chi_2()>0){
	    mu_id1 = muID1.Muon_likelihood();}
	  if (muID2.Chi_2()>0){
	    mu_id2 = muID2.Muon_likelihood();}

	  //	    	    cout << "muid1 " << mu_id1 << endl;
	  //	    	    cout << "muid2 " << mu_id2 << endl;
	  eid sel_e1(el1), sel_e2(el2);
	  double e_id1 = sel_e1.prob( 3, -1, 5 ); 
	  double e_id2 = sel_e2.prob( 3, -1, 5 ); 
	 
	  //	  cout << "EID " << e_id1 << " " << e_id2 << endl;

	  Pi1.GetImpactParameters(&dr1,&dz1,&nrp1,&nz1,&nax1,&nst1,0);
	  Pi2.GetImpactParameters(&dr2,&dz2,&nrp2,&nz2,&nax2,&nst2,0);

	  /*
	  nrp1 = trk_fit1.nhits(3);
	  nz1 = trk_fit1.nhits(4);
	  nrp2 = trk_fit2.nhits(3);
	  nz2 = trk_fit2.nhits(4);
	  */

	  /*	  cout << "SVD hits rp1 " << nrp1 << endl;
	  cout << "SVD hits z1 " << nz1 << endl;
	  cout << "SVD hits rp2 " << nrp2 << endl;
	  cout << "SVD hits z2 " << nz2 << endl;
	  */

	  kvertexfitter kv2pi;
	  int fitstatus_2pi = -1;

	  addTrack2fit(kv2pi, Part1);
	  addTrack2fit(kv2pi, Part2);
	  if (IpProfile::usable()) addTube2fit(kv2pi);
	  fitstatus_2pi = kv2pi.fit();
	  
	  /* cout << "kv2pi.fit() "  << kv2pi.fit() <<endl;
	     cout << "kv2pi.vertex() "  << kv2pi.vertex() <<endl;
	      cout << "kv2pi.errVertex() "  << kv2pi.errVertex() <<endl;
	      cout << "kv2pi.chisq() "  << kv2pi.chisq() <<endl;
	  */

	  float svd2pi =0;
	  svd2pi = 1000*nrp1+100*nz1+10*nrp2+nz2;
	  float fchisq2pi = kv2pi.chisq();
	  double prob2pi = -1.; 
	  if (IpProfile::usable()&&fchisq2pi>0.&&fitstatus_2pi==0) prob2pi=prob_(fchisq2pi, 3);
	  if (!IpProfile::usable()&&fchisq2pi>0.&&fitstatus_2pi==0) prob2pi=prob_(fchisq2pi, 1);
	  double pt_Pi1 = sqrt(Pi1.P4.px()*Pi1.P4.px() + Pi1.P4.py()*Pi1.P4.py());
	  double pt_Pi2 = sqrt(Pi2.P4.px()*Pi2.P4.px() + Pi2.P4.py()*Pi2.P4.py());

	  HepParticle P2PI(Pi1,Pi2);


              HepLorentzVector P_REC  = P_BEAM-Pi1.P4-Pi2.P4;
              HepLorentzVector P_REC1  = P_BEAM-Pi1.P4;
              HepLorentzVector P_REC2  = P_BEAM-Pi2.P4;
	      double m_rec  = P_REC.m();
	      double m_rec1  = P_REC1.m();
	      double m_rec2  = P_REC2.m();
	      HepLorentzVector P_REC_com = P_REC;

		  if (P_REC.m()>=9.8&&P_REC.m()<9.85){
		    if (my_npipi[0]<9){ my_npipi[0]=my_npipi[0]+1;}}
		  if (P_REC.m()>=9.85&&P_REC.m()<9.90){
		    if (my_npipi[1]<9){ my_npipi[1]=my_npipi[1]+1;}}
		  if (P_REC.m()>=9.90&&P_REC.m()<9.95){
		    if (my_npipi[2]<9){ my_npipi[2]=my_npipi[2]+1;}}
		  if (P_REC.m()>=9.95&&P_REC.m()<10.0){
		    if (my_npipi[3]<9){ my_npipi[3]=my_npipi[3]+1;}}
		  if (P_REC.m()>=10.0&&P_REC.m()<10.05){
		    if (my_npipi[4]<9){ my_npipi[4]=my_npipi[4]+1;}}
		  if (P_REC.m()>=10.05&&P_REC.m()<10.15){
		    if (my_npipi[5]<9){ my_npipi[5]=my_npipi[5]+1;}}
		  if (P_REC.m()>=10.15&&P_REC.m()<10.20){
		    if (my_npipi[6]<9){ my_npipi[6]=my_npipi[6]+1;}}
		  if (P_REC.m()>=10.20&&P_REC.m()<10.30){
		    if (my_npipi[7]<9){ my_npipi[7]=my_npipi[7]+1;}}
		  if (P_REC.m()>=10.30&&P_REC.m()<10.35){
		    if (my_npipi[8]<9){ my_npipi[8]=my_npipi[8]+1;}}
		  if (P_REC.m()>=10.35&&P_REC.m()<10.40){
		    if (my_npipi[9]<9){ my_npipi[9]=my_npipi[9]+1;}}


	      //	      cout << P_REC_com << endl;

	      P_REC_com.boost(-P_BEAM.boostVector());

	      /*	      cout << P_BEAM << endl;
	      cout << P_REC_com << endl;
	      cout << endl;
	      */
	      double hb_angle_com = P_REC_com.pz()/P_REC_com.rho();
	      double hb_angle = P_REC.pz()/P_REC.rho();

	      //	      cout << "hb_angle_com = " << hb_angle_com << endl;
	      //cout << "hb_angle = " << hb_angle << endl;

	      if (m_rec < 9.2 || m_rec > 11.0) continue;

	      if( fabs(dz1) < 5. && fabs(dz2) < 5. && fabs(dr1) < 1. && fabs(dr2) < 1. && //prob2pi >1.E-6 && svd2pi > 3332 &&
			 (Pi1.charged->charge()*Pi2.charged->charge()) < 0
		  &&  evt_R2<0.4
		  &&  evt_R2<0.35
		  &&  evt_R2<0.32
		  //&& (nax1+nst1)>9 && (nax2+nst2)>9
		  //		  && prob2pi > 1.E-4
		  )

		{ //if cuts

		  savecount++;
                      
		  *status = 1;
				 	         
		  //	       	  int my_ngam[10][10];
	       	  int my_ngam[10][8];
		  int gam_bini = 0;
		  int gam_binj = 0;
		  
		  //		  for (gam_bini=0;gam_bini<10;++gam_bini){
		  //for (gam_binj=0;gam_binj<10;++gam_binj){
		  for (gam_bini=0;gam_bini<8;++gam_bini){
		    for (gam_binj=0;gam_binj<10;++gam_binj){
		      my_ngam[gam_binj][gam_bini]=0;
		    }
		  }
		  //loop over it3

		  // double loop of i, j to fill 100 bins

		  // fill variable ngam0, ngam1

		  for(std::vector<HepParticle>::iterator it3 = myGamList.begin(); it3 != myGamList.end(); ++it3)
		    {
		      
		      HepParticle& G = *it3;
		      
		      HepLorentzVector P_GREC = P_REC-G.P4;
		      
		      HepLorentzVector P4_hbframe= G.P4;
		      
		      
		      double anglegwrthb = -10;
		      
		      
		      //	      P4_hbframe.boost((P2PI.P4).boostVector());
		      
		      P4_hbframe.boost(-P_REC.boostVector());
		      
		      //		      cout << "photon energy " << P4_hbframe.e() << endl;
		      //                      cout << "recoil energy " << P_REC.e() << endl;
		      double en=0.0;
		      double enstep = 0.0125;

		      //		      for ( int i=0; i < 10; i++) {
		      //for (int j=0; j < 10; j++) {		      
		      for ( int i=0; i < 10; i++) {
			for (int j=0; j < 8; j++) {		      
			  if (P4_hbframe.e()>=en&&P4_hbframe.e()<en+enstep){
			    if (my_ngam[i][j]<9){ 
			      ++my_ngam[i][j]; 
			    }}
			  en=en+enstep;
			}	      
		      }
		      
		      /*		      for ( int i=0; i < 10; i++) {
			for (int j=0; j < 10; j++) {		      
			  if( my_ngam[i][j] > 0 ) {
			  cout << "photon energy after " << P4_hbframe.e() << " " << i << " " << j << " " << my_ngam[i][j] << endl;
			  }
			}	      
		      }
		      */

		      /*
		      if (P4_hbframe.e()>=0.100&&P4_hbframe.e()<0.180){
			if (my_ngam[0]<9){ my_ngam[0]=my_ngam[0]+1;}}
		      if (P4_hbframe.e()>=0.050&&P4_hbframe.e()<0.230){
			if (my_ngam[1]<9){ my_ngam[1]=my_ngam[1]+1;}}
		      if (P4_hbframe.e()>=0.210&&P4_hbframe.e()<0.290){
			if (my_ngam[2]<9){ my_ngam[2]=my_ngam[2]+1;}}
		      if (P4_hbframe.e()>=0.160&&P4_hbframe.e()<0.340){
			if (my_ngam[3]<9){ my_ngam[3]=my_ngam[3]+1;}}
		      if (P4_hbframe.e()>=0.380&&P4_hbframe.e()<0.520){
			if (my_ngam[4]<9){ my_ngam[4]=my_ngam[4]+1;}}
		      if (P4_hbframe.e()>=0.330&&P4_hbframe.e()<0.570){
			if (my_ngam[5]<9){ my_ngam[5]=my_ngam[5]+1;}}
		      if (P4_hbframe.e()>=0.420&&P4_hbframe.e()<0.560){
			if (my_ngam[6]<9){ my_ngam[6]=my_ngam[6]+1;}}
		      if (P4_hbframe.e()>=0.370&&P4_hbframe.e()<0.610){
			if (my_ngam[7]<9){ my_ngam[7]=my_ngam[7]+1;}}
		      if (P4_hbframe.e()>=0.700&&P4_hbframe.e()<0.900){
			if (my_ngam[8]<9){ my_ngam[8]=my_ngam[8]+1;}}
		      if (P4_hbframe.e()>=0.900){
			if (my_ngam[9]<9){ my_ngam[9]=my_ngam[9]+1;}}
		      
		      */		      
		      
		      /*
			cout << "P_BEAM=" << P_BEAM <<endl;
			cout << "Pi1.P4=" << Pi1.P4 <<endl;
			cout << "Pi2.P4=" << Pi2.P4 <<endl;
			cout << "G.P4=" << G.P4 <<endl;
		      */
		      
		      
		      /*
			cout << "gammas=" << GamList.size() <<endl;
			cout << "tracks=" << TrkList.size() <<endl;
			cout << "pi0=" << pi0_count <<endl;
		      */
		      
		      
		      //	      if(G.P4.e() > 0.05 && GamList.size() < 20 && pi0_count < 10 && TrkList.size() < 20)
		      double m_grec = P_GREC.m();
		      
		      if (m_rec < 9.2 || m_rec > 11.0) continue;
		      
		      if( fabs(dz1) < 5. && fabs(dz2) < 5. && fabs(dr1) < 1. && fabs(dr2) < 1. && //prob2pi >1.E-6 && svd2pi > 3332 &&
			  (Pi1.charged->charge()*Pi2.charged->charge()) < 0
			  &&  evt_R2<0.4
			  //&& (nax1+nst1)>9 && (nax2+nst2)>9
			  //		  && prob2pi > 1.E-4
			  )
			
			{
			  
			  
			  savecount++;
			  
			  *status = 1;
			  
			}
		      
		    }	    

		  //fill ngam0, ngam1 ete here
		  int my_ngam_count0=0;
		  int my_ngam_count1=0;
		  int my_ngam_count2=0;
		  int my_ngam_count3=0;
		  int my_ngam_count4=0;
		  int my_ngam_count5=0;
		  int my_ngam_count6=0;
		  int my_ngam_count7=0;
		  int my_ngam_count8=0;
		  int my_ngam_count9=0;
		  
		  for (gam_bini=0;gam_bini<8;++gam_bini){
		    my_ngam_count0=my_ngam_count0+pow(10,gam_bini)*my_ngam[0][gam_bini];
		    my_ngam_count1=my_ngam_count1+pow(10,gam_bini)*my_ngam[1][gam_bini];
		    my_ngam_count2=my_ngam_count2+pow(10,gam_bini)*my_ngam[2][gam_bini];
		    my_ngam_count3=my_ngam_count3+pow(10,gam_bini)*my_ngam[3][gam_bini];
		    my_ngam_count4=my_ngam_count4+pow(10,gam_bini)*my_ngam[4][gam_bini];
		    my_ngam_count5=my_ngam_count5+pow(10,gam_bini)*my_ngam[5][gam_bini];
		    my_ngam_count6=my_ngam_count6+pow(10,gam_bini)*my_ngam[6][gam_bini];
		    my_ngam_count7=my_ngam_count7+pow(10,gam_bini)*my_ngam[7][gam_bini];
		    my_ngam_count8=my_ngam_count8+pow(10,gam_bini)*my_ngam[8][gam_bini];
		    my_ngam_count9=my_ngam_count9+pow(10,gam_bini)*my_ngam[9][gam_bini];
		  }
		  
		  int my_npipi_count=0;
		  
		  for (pipi_bin=0;pipi_bin<10;++pipi_bin){
		    my_npipi_count=my_npipi_count+pow(10,pipi_bin)*my_npipi[pipi_bin];
		  }
		  
		  m_tpl->column("ntr",     TrkList.size());
		  m_tpl->column("npi0",    pi0_count);
		  m_tpl->column("ngam",    myGamList.size());
		  m_tpl->column("npostr",  Npos);
		  m_tpl->column("nnegtr",  Nneg);
		  m_tpl->column("qpi1",    Pi1.charged->charge());
		  m_tpl->column("qpi2",    Pi2.charged->charge());
		  m_tpl->column("rp1",   nrp1);
              m_tpl->column("rz1",   nz1);
              m_tpl->column("rp2",   nrp2);
              m_tpl->column("rz2",   nz2);
              m_tpl->column("drpi1",   dr1);
              m_tpl->column("drpi2",   dr2);
              m_tpl->column("dzpi1",   dz1);
              m_tpl->column("dzpi2",   dz2);
              m_tpl->column("kprobpi1",kid.prob(Pi1.charged));              
	      m_tpl->column("kprobpi2",kid.prob(Pi2.charged));
	      m_tpl->column("plabpi1", Pi1.P4.rho());
	      m_tpl->column("plabpi2", Pi2.P4.rho());
	      /*
              m_tpl->column("pxpi1",   Pi1.P4.px());
              m_tpl->column("pypi1",   Pi1.P4.py());
              m_tpl->column("pzpi1",   Pi1.P4.pz());
              m_tpl->column("pxpi2",   Pi2.P4.px());
              m_tpl->column("pypi2",   Pi2.P4.py());
              m_tpl->column("pzpi2",   Pi2.P4.pz());
              m_tpl->column("pxg",     G.P4.px());
              m_tpl->column("pyg",     G.P4.py());
	      */
	      //	      m_tpl->column("pzg",    P4_hbframe.pz()); //G.P4_boost.pz());
	      //m_tpl->column("egboost", P4_hbframe.e()); // G.P4_boost.e());
	      //              m_tpl->column("eg",      G.P4.e());
              m_tpl->column("recm2pi", m_rec);
	      //              m_tpl->column("recmg2pi",m_grec);
              m_tpl->column("m2pi",    P2PI.P4.m());
              m_tpl->column("angle2pi",    hb_angle_com);
	      m_tpl->column("ebeam",   E_BEAM);
	      if (EvtMgr.count()>0) 
		{
		  m_tpl->column("expno",EvtMgr[0].ExpNo());
		  m_tpl->column("runno",EvtMgr[0].RunNo());
		  m_tpl->column("evtno",EvtMgr[0].EvtNo()&0x0FFFFFFF);
		}
	      m_tpl->column("ncls",(float)evt_Ncls);
	      m_tpl->column("nclsb",(float)Necl_barrel);
	      m_tpl->column("psum",evt_Psum);
	      m_tpl->column("esum",evt_Esum);
	      m_tpl->column("evis",evt_Evis);
	      m_tpl->column("evtpz",evt_Pz);
	      m_tpl->column("hjmass",evt_HeavyJetMass);
	      m_tpl->column("thrust",evt_Thrust);
	      m_tpl->column("evtr2",evt_R2);
	      m_tpl->column("chisq",fchisq);
	      m_tpl->column("chisqsv",fchisqsv);
	      m_tpl->column("chisq2pi",kv2pi.chisq());

	      //	      m_tpl->dumpData();

	      if (EvtMgr.count()>0) 
		{
		  m_new->column("expno",EvtMgr[0].ExpNo());
		  m_new->column("runno",EvtMgr[0].RunNo());
		  m_new->column("evtno",EvtMgr[0].EvtNo()&0x0FFFFFFF);
		}
	      m_new->column("ntr",     TrkList.size());
	      m_new->column("npi0",    pi0_count);
	      //	      m_new->column("ngam",    myGamList.size());
	      m_new->column("ngam0",    my_ngam_count0);
	      m_new->column("ngam1",    my_ngam_count1);
	      m_new->column("ngam2",    my_ngam_count2);
	      m_new->column("ngam3",    my_ngam_count3);
	      m_new->column("ngam4",    my_ngam_count4);
	      m_new->column("ngam5",    my_ngam_count5);
	      m_new->column("ngam6",    my_ngam_count6);
	      m_new->column("ngam7",    my_ngam_count7);
	      m_new->column("ngam8",    my_ngam_count8);
	      m_new->column("ngam9",    my_ngam_count9);

	      /*
	      cout << "Inside nutplemaker where m2mu = -1" << endl;
	      cout << "--------" << endl;
		  cout << "my_ngam_count0 = " << my_ngam_count0 << endl;
		  cout << "my_ngam_count1 = " << my_ngam_count1 << endl;
		  cout << "my_ngam_count2 = " << my_ngam_count2 << endl;
		  cout << "my_ngam_count3 = " << my_ngam_count3 << endl;
		  cout << "my_ngam_count4 = " << my_ngam_count4 << endl;
		  cout << "my_ngam_count5 = " << my_ngam_count5 << endl;
		  cout << "my_ngam_count6 = " << my_ngam_count6 << endl;
		  cout << "my_ngam_count7 = " << my_ngam_count7 << endl;
		  cout << "my_ngam_count8 = " << my_ngam_count8 << endl;
		  cout << "my_ngam_count9 = " << my_ngam_count9 << endl;
		  cout << "---------" << endl;
	      */
	      //              m_new->column("qpi1",    Pi1.charged->charge());
              //m_new->column("qpi2",    Pi2.charged->charge());
	      m_new->column("energy",E_BEAM);
	      //	      m_new->column("chisq2pi",kv2pi.chisq());
	      //m_new->column("probev",probev);
	      m_new->column("prob2pi",prob2pi);
	      //	      m_new->column("chi2ev",fchisqsv);
	      //m_new->column("chi22pi",fchisq2pi);
	      //	      m_new->column("svd2pi",svd2pi);
	      m_new->column("angle2pi",hb_angle_com);
              m_new->column("kprobpi1",kid.prob(Pi1.charged));              
	      m_new->column("kprobpi2",kid.prob(Pi2.charged));
              m_new->column("pprobpi1",prid.prob(Pi1.charged));              
	      m_new->column("pprobpi2",prid.prob(Pi2.charged));
	      m_new->column("eid1",e_id1);
	      m_new->column("eid2",e_id2);

	      m_new->column("muid1",mu_id1);
	      m_new->column("muid2",mu_id2);
	      m_new->column("evtr2",evt_R2);
              m_new->column("recm2pi", m_rec);
              m_new->column("m2pi",    P2PI.P4.m());
              m_new->column("pqpi1", (Pi1.charged->charge())*Pi1.P4.rho());
              m_new->column("pqpi2", (Pi2.charged->charge())*Pi2.P4.rho());


	      m_new->column("m2mu", -1);
	      //	      m_new->column("recmg2pi", m_grec);
	      //	      m_new->column("eg",  G.P4.e());
	      //m_new->column("egboost",  P4_hbframe.e());
	      //anglegwrthb = (P_REC.px()*P4_hbframe.px()+P_REC.py()*P4_hbframe.py()+P_REC.pz()*P4_hbframe.pz())/P_REC.rho()*P4_hbframe.rho();
	      //m_new->column("angleg",  anglegwrthb);
	      //m_new->column("anglelabg",G.P4.pz()/G.P4.rho());

	      m_new->column("rp1",   nrp1);
              m_new->column("rz1",   nz1);
              m_new->column("rp2",   nrp2);
              m_new->column("rz2",   nz2);
	      //              m_new->column("ax1",   nax1);
              //m_new->column("st1",   nst1);
              //m_new->column("ax2",   nax2);
              //m_new->column("st2",   nst2);
              m_new->column("dr1",   dr1);
              m_new->column("dr2",   dr2);
              m_new->column("dz1",   dz1);
              m_new->column("dz2",   dz2);


              m_new->column("recmpi1", m_rec1);
              m_new->column("recmpi2", m_rec2);

	      //m_new->column("npostr",  Npos);
	      //m_new->column("nnegtr",  Nneg);

	      //m_new->column("mchrg", (P_CHG-Pi1.P4-Pi2.P4).m());

	      m_new->column("npipi",    my_npipi_count);


	      m_new->dumpData();

	      if (m_rec1 < 10.55 && m_rec2 < 10.55) continue;
	      if (m_rec1 > 10.75 && m_rec2 > 10.75) continue;

	      
	      
	      if( //G.P4.e() > 0.2 && P4_hbframe.e() > 0.3 &&
		 fabs(dz1) < 5. && fabs(dz2) < 5. && fabs(dr1) < 1. && fabs(dr2) < 1 &&
		 fabs(dz1) < 2. && fabs(dz2) < 2. && fabs(dr1) < 0.3 && fabs(dr2) < 0.3 &&
		 (Pi1.charged->charge()*Pi2.charged->charge()) < 0
		 && evt_R2<0.4
		 && evt_R2<0.35
		 && evt_R2<0.3
		 && (nax1+nst1)>9 && (nax2+nst2)>9
		 //		 && prob2pi > 1.E-4
		 && prob2pi > 1.E-6
		 && e_id1< 0.9 && e_id2<0.9
		 && kid.prob(Pi1.charged)<0.9 && kid.prob(Pi2.charged)<0.9                           
		 && prid.prob(Pi1.charged)<0.9 && prid.prob(Pi2.charged)<0.9                           
 		  //				 && prob2pi >1.E-6  && svd2pi > 3332
				 // && P4_hbframe.e() > 0.3
		 )
		{

	      // add gamma loop here
	      for(std::vector<HepParticle>::iterator it3 = myGamList.begin(); it3 != myGamList.end(); ++it3)
		{
		  
		  HepParticle& G = *it3;
		  
		  HepLorentzVector P_GREC = P_REC-G.P4;
		  
		  HepLorentzVector P4_hbframe= G.P4;
		  
		  
		  double anglegwrthb = -10;
		  
		  double m_grec = P_GREC.m();
		  
		  //	      P4_hbframe.boost((P2PI.P4).boostVector());
		  
		  P4_hbframe.boost(-P_REC.boostVector());

		  if (m_grec < 9.2 || m_grec > 11.0) continue;
		  //		  if (m_rec < 9.8 || m_rec > 10.5) continue;
		  if (m_rec < 9.8 || m_rec > 10.35) continue;
		  if (P4_hbframe.e() < 0.1) continue;		  

		  m_new->column("m2mu", 0); //use this to index photon
		  m_new->column("recmg2pi", m_grec);
		  m_new->column("eg",  G.P4.e());
		  m_new->column("egboost",  P4_hbframe.e());
		  anglegwrthb = (P_REC.px()*P4_hbframe.px()+P_REC.py()*P4_hbframe.py()+P_REC.pz()*P4_hbframe.pz())/P_REC.rho()*P4_hbframe.rho();
		  m_new->column("angleg",  anglegwrthb);
		  m_new->column("anglelabg",G.P4.pz()/G.P4.rho());
		  
	      if (EvtMgr.count()>0) 
		{
		  m_new->column("expno",EvtMgr[0].ExpNo());
		  m_new->column("runno",EvtMgr[0].RunNo());
		  m_new->column("evtno",EvtMgr[0].EvtNo()&0x0FFFFFFF);
		}
	      m_new->column("ntr",     TrkList.size());
	      m_new->column("npi0",    pi0_count);
	      //	      m_new->column("ngam",    myGamList.size());
	      m_new->column("ngam0",    my_ngam_count0);
	      m_new->column("ngam1",    my_ngam_count1);
	      m_new->column("ngam2",    my_ngam_count2);
	      m_new->column("ngam3",    my_ngam_count3);
	      m_new->column("ngam4",    my_ngam_count4);
	      m_new->column("ngam5",    my_ngam_count5);
	      m_new->column("ngam6",    my_ngam_count6);
	      m_new->column("ngam7",    my_ngam_count7);
	      m_new->column("ngam8",    my_ngam_count8);
	      m_new->column("ngam9",    my_ngam_count9);

	      /*	      
	      cout << "Inside nutplemaker where m2mu = 0" << endl;
	      cout << "--------" << endl;
	      cout << "egboost: " << P4_hbframe.e() << endl;
	      cout << "recm2pi: " << m_rec << endl;
	      cout << "--------" << endl;
		  cout << "my_ngam_count0 = " << my_ngam_count0 << endl;
		  cout << "my_ngam_count1 = " << my_ngam_count1 << endl;
		  cout << "my_ngam_count2 = " << my_ngam_count2 << endl;
		  cout << "my_ngam_count3 = " << my_ngam_count3 << endl;
		  cout << "my_ngam_count4 = " << my_ngam_count4 << endl;
		  cout << "my_ngam_count5 = " << my_ngam_count5 << endl;
		  cout << "my_ngam_count6 = " << my_ngam_count6 << endl;
		  cout << "my_ngam_count7 = " << my_ngam_count7 << endl;
		  cout << "my_ngam_count8 = " << my_ngam_count8 << endl;
		  cout << "my_ngam_count9 = " << my_ngam_count9 << endl;
		  cout << "---------" << endl;

	      */
	      //              m_new->column("qpi1",    Pi1.charged->charge());
              //m_new->column("qpi2",    Pi2.charged->charge());
	      m_new->column("energy",E_BEAM);
	      //	      m_new->column("chisq2pi",kv2pi.chisq());
	      //	      m_new->column("probev",probev);
	      m_new->column("prob2pi",prob2pi);
	      //	      m_new->column("chi2ev",fchisqsv);
	      //m_new->column("chi22pi",fchisq2pi);
	      //m_new->column("svd2pi",svd2pi);
	      m_new->column("angle2pi",hb_angle_com);
              m_new->column("kprobpi1",kid.prob(Pi1.charged));              
	      m_new->column("kprobpi2",kid.prob(Pi2.charged));
              m_new->column("pprobpi1",prid.prob(Pi1.charged));              
	      m_new->column("pprobpi2",prid.prob(Pi2.charged));
	      m_new->column("eid1",e_id1);
	      m_new->column("eid2",e_id2);

	      m_new->column("muid1",mu_id1);
	      m_new->column("muid2",mu_id2);
	      m_new->column("evtr2",evt_R2);
              m_new->column("recm2pi", m_rec);
              m_new->column("m2pi",    P2PI.P4.m());
              m_new->column("pqpi1", (Pi1.charged->charge())*Pi1.P4.rho());
              m_new->column("pqpi2", (Pi2.charged->charge())*Pi2.P4.rho());


	      //	      m_new->column("m2mu", -1);
	      //	      m_new->column("recmg2pi", m_grec);
	      //	      m_new->column("eg",  G.P4.e());
	      //m_new->column("egboost",  P4_hbframe.e());
	      //anglegwrthb = (P_REC.px()*P4_hbframe.px()+P_REC.py()*P4_hbframe.py()+P_REC.pz()*P4_hbframe.pz())/P_REC.rho()*P4_hbframe.rho();
	      //m_new->column("angleg",  anglegwrthb);
	      //m_new->column("anglelabg",G.P4.pz()/G.P4.rho());

	      m_new->column("rp1",   nrp1);
              m_new->column("rz1",   nz1);
              m_new->column("rp2",   nrp2);
              m_new->column("rz2",   nz2);
	      //              m_new->column("ax1",   nax1);
              //m_new->column("st1",   nst1);
              //m_new->column("ax2",   nax2);
              //m_new->column("st2",   nst2);
              m_new->column("dr1",   dr1);
              m_new->column("dr2",   dr2);
              m_new->column("dz1",   dz1);
              m_new->column("dz2",   dz2);


              m_new->column("recmpi1", m_rec1);
              m_new->column("recmpi2", m_rec2);

	      //m_new->column("npostr",  Npos);
	      //m_new->column("nnegtr",  Nneg);

	      //m_new->column("mchrg", (P_CHG-Pi1.P4-Pi2.P4).m());

	      m_new->column("npipi",    my_npipi_count);



		  
	      m_new->dumpData();
		}
		}		
		}
	  nit2++; eit2++;
	    }
	  nit1++; eit1++;
	}
  eventcount++;
  if (eventcount % 100==0) printf("Event selection - %d/%d\n",savecount,eventcount); 
}

/*
 .......................................................................
 . FWJET2 - Subroutine to compute event sphericity 
 . 
 . Inputs    : NTRK -- number of tracks in momentum array 
 .           : PTRK -- array dimensioned (3,NTRK) containing 3-momenta 
 .           :         of tracks in the event 
 .           : EBEAM -- beam energy of the event 
 . Outputs   : H0 -- 0th moment of the event 
 .           : R1 -- H1/H0 of the event 
 .           : R2 -- H2/H0 of the event 
 .           : R3 -- H3/H0 of the event 
 .           : R4 -- H4/H0 of the event 
 . 
 . Calls     : LEGN(X) -- INTERNAL functions for Legendre polynomials 
 . Called    : Anyone who wants the Fox/Wolfram moments of the event 
 . Author    : F. Morrow  09/02/89  14.47.58 
 .             Transfer to C by kfjack 2/9,1999
 . 
 . Detailed description 
 .     FWJET2 determines the Fox/Wolfram moments of the event, up to 
 .     order 5. 
 .......................................................................
*/
void hb::fwjet2(int NTrk, float *PTrk, float EBeam, 
		     float *H0, float *R1, float *R2, float *R3, float *R4)
{
  int i, j;
  float sss ,cxang, ppari, pparj, delh0;
  float h, h1, h2, h3, h4;

  if (NTrk >= 2) {
    
    if (EBeam < 4.0 || EBeam > 6.5) EBeam = 5.29;
    sss = EBeam*EBeam*4;
    
    h=h1=h2=h3=h4=0;
    
    for (i=0;i<NTrk;i++) {
      ppari = sqrt(PTrk[i*3]*PTrk[i*3]+
		                PTrk[i*3+1]*PTrk[i*3+1]+
		   PTrk[i*3+2]*PTrk[i*3+2]);
      
      for (j=i;j<NTrk;j++) {
	
	if (i!=j) {
	  pparj = sqrt(PTrk[j*3]*PTrk[j*3]+
		                    PTrk[j*3+1]*PTrk[j*3+1]+
		       PTrk[j*3+2]*PTrk[j*3+2]);

	  cxang = (PTrk[i*3]*PTrk[j*3]+ 
		    PTrk[i*3+1]*PTrk[j*3+1]+
		   PTrk[i*3+2]*PTrk[j*3+2])/(ppari * pparj);
	  
	  delh0 = 2*ppari*pparj;
	  
	  h += delh0;
	  h1 += delh0 * cxang;
	  h2 += delh0 * LGNDR2(cxang);
	  h3 += delh0 * LGNDR3(cxang);
	  h4 += delh0 * LGNDR4(cxang);
	} else {
	  delh0 = ppari*ppari;
	  h += delh0;
	  h1 += delh0;
	  h2 += delh0;
	  h3 += delh0;
	  h4 += delh0;
	}
      }
    }
  }

  if (h > 1e-7f) {
    *H0 = h / sss;
    *R1 = h1 / h;
    *R2 = h2 / h;
    *R3 = h3 / h;
    *R4 = h4 / h;
  }else {
    *H0=*R1=*R2=*R3=*R4=-1.0;
  }
}


#if defined(BELLE_NAMESPACE)
} // namespace Belle
#endif

